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1. Introduction 


Several times in science and engineering we need to evaluate rather exotic functions -- elliptic 
integrals, Bessel functions, and probability density functions, for example. Here’s where you can get 
your hands on some LISP code that will evaluate some of these functions. The code runs in both 
MACLISP and LISPM IJSP. 

The functions in the library are scattered in several different files. The descriptions below will 
describe the function name, how it should be used, and the file where that function can be obtained. For 
example, the function FACTOR is in PS:<GLR.FUNCT>MOD.LSP. Before it can be used, the code 
must first be loaded by doing one of: 


(LOAD "PS:<GLR.FUNCT>M0D.FASL") 

(LOAD "OZ:PS:<GLR.FUNCT>M0D.QFASL") 
(LOAD "OZ:PS:<GLR.FUNCT>M0D.QBIN") 

(LOAD "OZ:PS:<GLR.FUNCT>M0D.BIN" ) 


For MACLISP 
For MIT LISPM 
For Symbolics LM2 
For Symbolics 3600 


After the file has been loaded, just call the function with the right arguments: 


(FACTOR 12345678987654321) ==> (3 3 3 3 37 37 333667 333667) 

If you have to compile one of the library files or you want to read the source version of a library 
file into your LISP, then you will have to load a special macro called IMPORT-FILE that is defined in 
the file PS:<GLR.LISP>IMPORT.LSP. You may want to use the IMPORT-FILE macro yourself. Its 
format is: 


(IMPORT-FILE "OZ:PS:<GLR.FUNCT>M0D.LSP") 

The argument is die source filename widi a LISPM style hostname (the MACLISP version of die macro 
knows to ignore the host). The behavior of die macro depends on whether it is being compiled, 
evaluated, or loaded. If the IMPORT-FILE form is evaluated or loaded (EVAL-WHEN (EVAL LOAD) 
...) and die imported file has not already been loaded, it will load a binary file if it can find one or the 
source file if it cannot. At compile time (EVAL-WHEN (COMPILE)...), die IMPORT-FILE form turns 
into 


(PROGN ’COMPILE <first form in source file> <load code>) 

By convention, die first form in die imported source file should be a DECLARE of all die external 
functions and variables that die imported file defines. Including the declarations allows some compile 
time error checking and eliminates spurious warnings about undefined functions. The <loadcode> 
conspires to load die binary or source file when when the file we are currently compiling gets loaded. 
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2. Number Theory and Combinatorial Functions 

(FACTORIAL N) 

The number of ways to permute N objects. 


PS:<GLR.FUNCT>COMBIN.LSP 


(CHOOSE N M) PS:<GLR.FUNCT>COMBIN.LSP 

The number of ways to place M distinguishable objects on N distinguishable plates. 

(BELL-NUMBER N) PS:<GLR.FUNCT>COMBIN.LSP 

The number of ways to place N distinguishable objects on N indistinguishable plates. 

(CATALAN-NUMBER N) ps:<glr.funct>combin.lsp 

The number of ways to fully parendiesize a string of n symbols. 

(FIB N) PS: <GLR. FUNCT>FIB. LSP 

The Nth (FIXNUM) Fibonacci number. FIB(0) = 0, FIB(1) = 1,.... Uses a log(N) algorithm. 

(FACTOR N) PS : <GLR . FUNCT>M0D. LSP 

Find the prime factors of the integer N (FIXNUM or BIGNUM). Returns a sorted list of the 

factors. 


(PRIME-TEST N Optional (TRIALS 50.)) PS:<GLR.FUNCT>M0D.LSP 

Tests the primality of N using a probabilistic algorithm (see <Solovay>). Returns NIL, PRIME, 
or PROBABLY-PRIME (P{crror} = 2' TR,ALS ). TRIALS must be FIXNUM. 

(SMALLER-PRIME N) ps:<glr.funct>mod.lsp 

Finds the first prime that is smaller than N. Repeatedly calls PRIME-TEST. 


(JACOBI-SYMBOL P Q) 

Computes the Jacobi Symbol of P and Q. 


PS:<GLR.FUNCT>M0D.LSP 


(TOTIENT N) PS : <GLR. FUNCT>M0D. LSP 

Computes Euler’s totient function of N. Calls FACTOR as a subroutine. 


3. Discrete Fourier Transform 

These functions are used to compute discrete fourier transforms <Oppcnhcimer>. The inputs 
are two FLONUM arrays that represent the real and imaginary parts of the input. The algorithms require 
input arrays for the transform to be bit reversed; there are functions for doing die reversal. The functions 
also require some initialization. 'Hie length of the DPI’ must be a power of 2. 
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3 OISCRFTF FOURIFR TRANSFORM 


X[k] = X CX P(' j 2 77 k m / N) 
x[k] = (1/N) 2 N;1 q cxp( j 2 77 k m / N) 


(DFT-INIT LOGN) PS:<GLR.FUNCT>DFT.LSP 

DFT-INIT returns a structure containing some tables that are needed by the DFT algorithm. 
This structure should be saved away and given to the DFT routine each time it is called. The structure 
need only be computed once for each size DFT. The length of the DPT is 2 ^GN 

(DFT-FORWARD X-REAL X-IMAG TABLES) PS:<GLR.FUNCT>DFT.LSP 

This function docs a discrete Fourier transform. The results arc written back into X-REAL and 
X-IMAG (the previous contents are lost) and are in sequential order (ie -- not bit reversed). TABLES is 
as returned by DPT-INIT. 

(DFT-REVERSE X-REAL X-IMAG TABLES) PS:<GLR.FUNCT>DF T.LSP 

This is just like DFT-FORWARD, but it does the inverse transform. The input X arrays should 
be in sequential (not bit-reversed) order. The resulting arrays are also in sequential order. 

(DFT-REVERSE-ARRAY ARRAY TABLES) PS:<GLR.FUNCT>DFT.LSP 

This bit reverses the first 2 ^GN elements of ARRAY. 

(DFT-610 X-REAL X-IMAG TABLES) PS:<GLR.FUNCT>DFT.LSP 

This function docs a forward transform. X-REAL and X-IMAG are bit-reversed input arrays 
(to reverse them, use DFT-REVERSP>ARRAY). The results are returned in X-REAL and X-IMAG in 
sequential order. 

(DFT-618 X-REAL X-IMAG TABLES) PS:<GLR.FUNCT>DFT.LSP 

The reverse transform. Input arrays are sequential order, output arrays arc bit-reversed. Use 
DFT-REVERSE-ARRAY to put them in sequential order. 

4. Trigonometric and other Functions 

These approximations come from <Abramowitz>, <DEC>, and <Hastings>. 

(EXP10 X) PS:<GLR.FUNCT>EXTFCN.LSP 

Computes lOtX. 
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(LOGIO X) PS:<GLR.FUNCT>EXTFCN.LSP 

Log bclSC 10 of X. 

(TAN X) 

(SEC X) 

(CSC X) 

(COT X) 

(ASIN X) 

(ACOS X) 

(ASEC X) 

(ACSC X) 

(ACOT Y X) PS:<GLR.FUNCT>HYPER.LSP 

Various trigonometric functions. 

(ERROR-FUNCTION X) 

Error function. eps< 1.5E-7. 

(2/sqrt(77)) /q X exp(-t^) dt 

(BESSEL-I N X) 

Modified Bessel function for integer order N. eps < 1.6E-7 

(BESSEL-J N X) 

Bessel function for integer order N. eps < 5E-8. 

(GAMMA-FUNCTION A) ps:<glr.funct>gamma.lsp 

Gamma function. 

T(a) = /q°° t 3 ’ 1 c _t dt 

(GAMMA-FUNCTION-INCOMPLETE A X) ps : <glr.funct>gamma.lsp 

Incomplete Gamma function. 0 <= x < = inf, but breaks if X is large because of roundoff error. 

T(a, x) - Jq x t a '* e _t dt 

(BETA-FUNCTION A B) ps:<glr.funct>gamma.lsp 

Beta function. A and B arc FLONUM. 

fMa. b) = /q 1 t a_1 (l-t) b_1 dt, a>0, b>0 

(BETA-FUNCTION-INCOMPLETE A B X) ps:<glr.funct>gamma.lsp 

Incomplete Beta function. A and B are FLONUM, but one of A or B must be an integer (eg, 

/?(a, b, x) = / 0 X t a_1 (l-t) b_1 dt, (Ka, (Kb, (X = x< = 1 


PS:<GLR.FUNCT>EXTFCN.LSP 


PS:<GLR.FUNCT>BESSEL.LSP 


PS:<GLR.FUNCT>BESSEL.LSP 


3.0). 
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4 TRIGONOMETRIC AND OTHER EUNCTIONS 


(ELLIPTIC-INTEGRAL-K M) PS:<6LR.FUNCT>ELLIP.LSP 

Elliptic integral of the first kind, k(M). 0 <= M < 1. eps = 2.0E-8. 

f 0 7T/2 (1 - m sin 2 (0))' 0 - 5 dO 

(ELLIPTIC-INTEGRAL-KC M) PS:<GLR.FUNCT>ELLIP.LSP 

Complementary elliptic integral of the first kind, k’(M). k’(M) = k(l-M). 

(ELLIPTIC-INTEGRAL-E M) PS:<GLR.FUNCT>ELLIP.LSP 

Elliptic integral of the second kind, c(M). eps < 2.0E-8. 

Jq 7772 sqrt(l - m sin 2 (0)) dO 

(ELLIPTIC-INTEGRAL-EC M) PS:<GLR.FUNCT>ELLIP.LSP 

Complementary elliptic integral of die second kind, e’(M). e’(M) = e(l-M). 

(ELLIPTIC-SINE u M) PS:<GLR.FUNCT>ELLIP.LSP 

Elliptic sine function (SN(u,M)) 

(ELLIPTIC-COSINE u M)) PS:<GLR.FUNCT>ELLIP.LSP 

Elliptic cosine function (CN(u,M)) 

5. Linear Regression 

Say we are trying to fit some (x, y) data to a function that looks like: 

Y’(x) = aO* 1 + al*gl(x) + a2*g2(x) + ... + an*gn(x) 

where the a[i] arc constants that we are trying to determine, the g[i] are linearly independent functions 
that we specify, and Y’(x) is the value of the fitted equation. Such a fit is called a linear regression and 
here arc some functions that will do it. 

(FIT FCTN X Y N Optional (MODE 0)) PS:<GLR.FUNCT>FIT.LSP 

Let X be a one dimensional array of x values (the x values could themselves be vectors). Y is a 
one dimensional array of FLONUM values for the corresponding X input. N is the number of terms in 
die equation we are fitting (ie, n in die equation above). FCTN is a function of N + l arguments that 
evaluates the equation above. That is FCTN looks like: 

(LAMBDA (X AO A1 A2 A3 ... AN) 

... ) 

and calculates die above equation. FIT returns a list of die coefficients Ai as the CDR of its value. 

FIT actually uses die procedure REGRESS in PS:<GLR.FUNCT>REGRES.LSP, which is 
modeled after <Bcvington>. If you want to do something a little more complicated, then ask me about 
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that procedure. 


6. Functional Minimization 

(FMFP FUNCT N X G EST EPS LIMIT) PS:<GLR.FUNCT>FMFP.LSP 

Flctcher-Powcll's functional minimization procedure (adapted from <Kuestcr> and 
<Lucnbcrgcr>). This finds a minimum of a multivariate, unconstrained, nonlinear function. That is, find 
the vector X such that F(X) is a local minimum. The arguments to this function are: 

N number of independent variables 

X vector[N] of initial variable values 

contains the result vector on exit 
G vector[N] in which to store the gradient 

EST estimate of minimum value of objective function 
EPS test value representing the expected absolute error in movement 
LIMIT maximum number of iterations 

FUNCT user supplied objective function that computes F and the gradient 
(FUNCT N X G) returns F(X) and 
also fills in the vector[N] G with the gradient 

(MARQUARDT N K X Y Z FUNC DERIV B BMIN BMAX BV) 

PS:<GLR.FUNCT>MARQ.LSP 

Marquardfs parameter fitting procedure (adapted from <Kuester>). Given N data points (X[i], 
Y[i]) and a function F with K parameters B[j], find the K parameters B[j] that best fit the data. The fitted 
values are returned in array Z. The arguments are: 

N -- number of data points (FIXNUM) 

K -- number of unknowns (FIXNUM) 

B -- vector[K] of (FLONUM) unknowns 

BMIN -- vector[K] of (FLONUM) minimum values of B[j] 

BMAX -- vector[K] of (FLONUM) maximum values of B[j] 

X -- vector[N] independent variable data points 
(this vector might be a vector of vectors) 

Y -- vector[N] of (FLONUM) dependent variable 

Z -- vector[N] of (FLONUM) computed values of dependent variable 

BV -- vector[K] of (FLONUM) codes 

1.0 -> numerical derivatives 
0.0 -> do not change this unknown 
FUNC -- (FUNC X K B N Z ZOFF) computes the function F (using the 
K parameters in vector B) for each of the X[0] ... X[N-1] 
and puts the respective results into Z[ZOFF+0] ... Z[Z0FF+N-1] 

DERIV-- (not used) 
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7 1IYP1IRB0LIC & INVFRSH HYPERBOLIC FUNCTIONS 


7. Hyperbolic & Inverse Hyperbolic Functions 

This functions use the Common Lisp <Stccle> names, but are only defined for real arguments. 
They call EXP to compute their results. 

(COSH X) 

(SINH X) 

(TANH X) 

(COTH X) 

(SECH X) 

(CSCH X) PS:<GLR.FUNCT>HYPER.LSP 

Hyperbolic functions. FLONUM argument and FLONUM result. 

(ACOSH X) 

(ASINH X) 

(ATANH X) 

(ACOTH X) 

(ASECH X) 

(ACSCH X) PS:<GLR.FUNCT>HYPER.LSP 

Inverse hyperbolic functions. FLONUM argument and FLONUM result 

8. Numerical Integration 

(INTEGRATE-TRAPEZOIDAL F XO XI N) PS:<GLR.FUNCT>INTEGR.LSP 

Uses the trapezoidal rule to integrate the function F (of one FLONUM argument) from XO to 
XI with N (FIXNUM) iterations. 

(INTEGRATE-SIMPSON F XO XI N) PS:<GLR.FUNCT>INTEGR.LSP 

Uses Simpson’s rule to integrate the function F (of one FLONUM argument) from XO to XI 
with N (FIXNUM) iterations. 

(INTEGRATE-EULER F XO YO H XI) PS:<GLR.FUNCT>RUNGE.LSP 

Uses the forward Euler method to integrate F(x, y) from position (XO, YO) to a position (XI, Yl) 
using a step size of H. Yl is the value of INTEGRATE-EULER. F is a function of X and Y (FLONUM) 
and should return DY/DX for the given position. 

(INTEGRATE-RUNGE-KUTTA F XO YO H XI) PS:<GLR.FUNCT>RUNGE.LSP 

This is just like EULER except it uses a fourth order RUNGE-KUITA method instead of 
Euler’s method. 
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(INTEGRATE-MULTI-STEP-H F XO YO H XI 

&optional (CD 1.0E-8) (CH 1.0E-5)) PS:<GLR.FUNCT>RUNGE.LSP 

Uses a multi-step predictor-corrector method to integrate DRRIV (sec <Mamming> page 407). 
Integrates from XO to XI using an initial step size of H. When die routine estimates that the error is 
below CD, the step size is doubled; when the error is above CH, the step size is halved. For 
reasonableness, CD « CH. 

9. Root Finding Functions 

(BISECTION-ROOT-SEARCH F XO XI) PS:<GLR.FUNCT>R00T.LSP 

Uses bisection search to find a zero of F(X) between XO and XI. 

(FALSE-POSITION-SEARCH F XO XI EPS) PS:<GLR.FUNCT>R00T.LSP 

Uses false position search to find a zero of F(X) with initial guesses XO and XI. 

(FALSE-POSITION-CONVERGE F XO XI EPS) PS:<GLR.FUNCT>ROOT.LSP 

Uses false position convergence to find an X = F(X) with initial guesses XO and XI. 

(CONVERGE F XO EPS) PS:<GLR.FUNCT>ROOT.LSP 

Find an X = F(X) with initial guess of XO. Uses Wegstcin’s method. 

10. Probability and Statistics 

The following functions are useful in probability and statistics. There are functions for 
computing probability density functions p(x), computing cumulative distributions P(x), and generating 
random numbers from a particular distribution. Most of die approximations come from <Abramowitz>. 
Most of the generators use die ratio of uniform deviates mediod <Kindcrman>. 

(UNIFORM-DENSITY X) 

(UNIFORM-CUMULATIVE X) 

(UNIFORM-RANDOM-NUMBER) PS:<GLR.FUNCT>STATIS.LSP 

Probability functions for die uniform distribution. 

p(x) = 1 0 <= x < — 1 

(NORMAL-DENSITY X) 

(NORMAL-CUMULATIVE X) PS:<GLR.FUNCT>STATIS.LSP 

Probability functions for normal (Gaussian) distribution. 

p(x) = (I/sqrt(2 7r))cxp (-x^/2) 
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10 PROBABILITY AND STATISTICS 


(NORMAL-TAIL ALPHA) PS:<GLR.FUNCT>STATIS.LSP 

Given a probability ALPHA, find the X such that ALPHA = l-P(X), where P(X) is the 
cumulative distribution. T his function provides a one-sided tail. 

(NORMAL-RANDOM-NUMBER) PS:<GLR.FUNCT>STATIS.LSP 

Generates a random number from the unit normal distribution. 

(EXPONENTIAL-DENSITY X LAMBDA) 

(EXPONENTIAL-CUMULATIVE X LAMBDA) 

(EXPONENTIAL-RANDOM-NUMBER LAMBDA) ps : <glr.funct>statis.lsp 

Exponential probability functions. 

p(x, X) = X exp(- X x), X > 0 

(POISSON-DENSITY N TIME) 

(POISSON-CUMULATIVE N TIME) 

(POISSON-RANDOM-NUMBER TIME) ps : <glr.funct>statis.lsp 

Poisson probability functions. Probability of exactly N (FIXNUM) arrivals within TIME 

(FLONUM) for a Poisson process with an average arrival rate (X) of 1. 

p(n, m) = (m n exp(-m))/n!, m = X t 

(CHI-SQUARE-DENSITY X N) 

(CHI-SQUARE-CUMULATIVE X N) 

(CHI-SQUARE-RANDOM-NUMBER N) ps:<glr.funct>statis.lsp 

The CHI-SQUARE distribution with N (FIXNUM) degrees of freedom. 

(T-DENSITY X N) 

(T-CUMULATIVE X N) 

(T-RANDOM-NUMBER N) ps:<glr.funct>statis.lsp 

Student’s T distribution with N (FIXNUM) degrees of freedom. N must be even in 

T-CUMULATIVE. 

(T-TWO-SIDED X N) ps:<glr.funct>statis.lsp 

For the T distribution, the probability that T falls within -X to +X. N must be an even 

FIXNUM. 

(F-DENSITY X M N) 

(F-CUMULATIVE X M N) 

(F-RANDOM-NUMBER M N) PS:<GLR.FUNCT>STATIS.LSP 

Snedccor’s F distribution with M and N (FIXNUM) degrees of freedom. 
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(GAMMA-DENSITY X A) 

(GAMMA-CUMULATIVE X A) 

(GAMMA-RANDOM-NUMBER A) ps:<glr.funct>statis.lsp 

Gamma probability functions with parameter A (FLONUM). The random number generator 
must have A > 1.0. 

p(x, a) = (l/r(a)) x 3 ' 1 e' x , a>0 


(BETA-DENSITY X A B) 

(BETA-CUMULATIVE X A B) 

(BETA-RANDOM-NUMBER A B) ps:<glr.funct>statis.lsp 

Beta probability functions with parameters A and B (FLONUM). The random number 
generator is slow, so keep A and B small (say below 10). 

p(x, a, b) = (l//?(a,b)) x a_1 (l-x) b_1 , a>0, b>0 


(CAUCHY-DENSITY X) 

(CAUCHY-CUMULATIVE X) 
(CAUCHY-RANDOM-NUMBER) 

Probability functions for the Cauchy distribution. 

p(x) = 1 / (7T (1 + X 2 )) 


PS:<GLR.FUNCT>STATIS.LSP 
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